Tutorial

Published

January 4, 2026

Introduction

This online tutorial walks you through each step of meta-analysis and location-scale meta-regression using three R packages: metafor, brms, and glmmTMB. With one real dataset, we show how to analyse mean effects, within-population variability, and context-dependent heterogeneity, from data preparation through model interpretation.

Detailed methodological explanations and additional examples are beyond the scope of this tutorial - readers interested in further details may find the following website and relevant papers helpful:

  1. Meta-analysis
    Online tutorial: https://itchyshin.github.io/Meta-analysis_tutorial/
    Related paper:
    Nakagawa S, Yang Y, Macartney EL, Spake R, Lagisz M. (2023). Quantitative evidence synthesis: a practical guide on meta-analysis, meta-regression, and publication bias tests for environmental sciences. Environmental Evidence. 12:8. https://doi.org/10.1186/s13750-023-00301-6

  2. Location–scale meta-analysis models
    Online tutorial: https://itchyshin.github.io/location-scale_meta-analysis/
    Related paper:
    Nakagawa S, Mizuno A, Morrison K, Ricolfi L, Williams C, Drobniak SM, Lagisz M, Yang Y. (2025). Location‐Scale Meta‐Analysis and Meta‐Regression as a Tool to Capture Large‐Scale Changes in Biological and Methodological Heterogeneity: A Spotlight on Heteroscedasticity. Global Change Biology. 31:e70204. https://doi.org/10.1111/gcb.70204

  3. Location–scale models
    Online tutorial: https://ayumi-495.github.io/Eco_location-scale_model/
    Related paper:
    Nakagawa S, Ortega S, Gazzea E, Lagisz M, Lenz A, Lundgren E, Mizuno A. (2025). Location–scale models in ecology and evolution: Heteroscedasticity in continuous, count and proportion data. Methods in Ecology and Evolution. https://doi.org/10.1111/2041-210x.70203

Note

Outputs for the scale component in scale parts are reported on the log scale. Guidance on how to interpret these parameters, including back-transformation and interpretation is provided in the third online tutorial.

Setup

Load packages

# needed packages
pacman::p_load("brms", "metafor", "tidyverse", 
               "tidybayes", "orchaRd", "here")

To do meta-analysis with glmmTMB, we need to install a development version of the package from GitHub:

# build package from branch: 
# probably remove this
remove.packages("glmmTMB")  ## if you have already installed glmmTMB

remotes::install_github("coraliewilliams/glmmTMB", ref = "equalto_covstruc", 
                        subdir = "glmmTMB", force = TRUE) 
# library(glmmTMB)
glmmTMB:::.valid_covstruct ## check the available covstruc (you should see equalto)

Data preparation

We use the dataset from Ponisio et al. (2015). First, we load the dataset and prepare it for the analyses.

datfull <- read.csv(here::here("data", "ponisio2014dataset.csv"), header = TRUE)
# dim(datfull)
# [1] 1071   72 -> 1071 effect sizes and 72 columns
# head(datfull) # check column names

# create ID using Author + Year + Journal - gives unique identifier for each study
datfull$study_ID <- as.numeric(as.factor(paste(datfull$Author, datfull$Year, datfull$Journal, sep = "_")))
datfull$effect_size_ID <- 1:nrow(datfull)

# filter to the analysis subset
# Org.fertilizer.type = either plant or animal and also Crop.type = cereals
dat <- datfull %>%
  filter(Org.fertilizer.type %in% c("plant", "animal")) %>%
  filter(Crop.type == "cereals") %>%
  droplevels()

# dim(dat) # 318 effect sizes

Effect size calculation

Original study calculated and used log response ratio (lnRR) and its variance - we need to get log variance ratio (lnVR) and log coefficient of variation ratio (lnCVR), respectively. metafor has a function escalc to calculate various effect sizes including lnVR and lnCVR.

dat <- escalc(measure = "VR",
                      m1i = Mean.Org,
                      sd1i = SD.Org,
                      n1i = N.Org,
                      m2i = Mean.conv,
                      sd2i = SD.conv,
                      n2i = N.Conv,
                      data = dat,
                      var.names = c("lnVR", "var.lnVR")
              )

dat <- escalc(measure = "CVR",
              m1i = Mean.Org,
              sd1i = SD.Org,
              n1i = N.Org,
              m2i = Mean.conv,
              sd2i = SD.conv,
              n2i = N.Conv,
              data = dat,
              var.names = c("lnCVR", "var.lnCVR")
              )

Model 1: Meta-analysis of mean

Okay, let’s start with a standard meta-analysis. We fit multilevel meta-analysis and meta-regression models to estimate the overall mean effect size (lnRR). Effect sizes were modelled as the response variable, with random intercepts specified for study identity and effect-size identity to account for the hierarchical structure of the data and non-independence among observations.

metafor

When we want to run the multilevel meta-analysis using metafor, we can use the rma.mv function. Here, we specify lnRR as the response variable and var.lnRR as the sampling variance. Random effects are specified using the random argument, where we include random intercepts for both study_ID and effect_size_ID.

mod_lnRR <- rma.mv(lnRR,
                    var.lnRR,
                    random = list(~ 1 | study_ID,
                                  ~ 1 | effect_size_ID),
                    data = dat,
                    test = "t",
                    method = "REML"
                   )
summary(mod_lnRR)
# Multivariate Meta-Analysis Model (k = 318; method: REML)
# 
#    logLik   Deviance        AIC        BIC       AICc   
# -109.8272   219.6545   225.6545   236.9312   225.7311   
# 
# Variance Components:
# 
#             estim    sqrt  nlvls  fixed          factor 
# sigma^2.1  0.0736  0.2712     36     no        study_ID 
# sigma^2.2  0.0830  0.2881    318     no  effect_size_ID 
# 
# Test for Heterogeneity:
# Q(df = 317) = 15841.3280, p-val < .0001
# 
# Model Results:
# 
# estimate      se     tval   df    pval    ci.lb    ci.ub      
#  -0.3204  0.0534  -5.9951  317  <.0001  -0.4255  -0.2152  *** 
# 
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

i2_ml(mod_lnRR)
         # I2_Total       I2_study_ID I2_effect_size_ID 
         # 98.11404          46.09489          52.01915 

The estimated average effect size was negative: \(\beta_0\) = -0.32, 95% confidence interval, CI = [-0.43, -0.22]), indicating that organic farming reduced crop yield by approximately 27% compared to conventional farming (exp(-0.32) = 0.73). There was also high heterogeneity among effect sizes (\(I^2_{\text{Total}}\) = 98.11%), with contribution s from both the study level and the effect size level (\(I^2_{\text{between}}\) = 46.09, \(I^2_{\text{within}}\) = 52.02). This indicates that effect sizes vary widely across studies as well as effect sizes…

Next, we fit a simple meta-regression to examine the type of fertiliser (animal vs plant) explains variation in mean effect sizes.

mod_lnRR_mr <- rma.mv(lnRR,
                    var.lnRR,
                    mods = ~ 1 + Org.fertilizer.type,
                    random = list(~ 1 | study_ID,
                                  ~ 1 | effect_size_ID),
                    data = dat,
                    test = "t",
                    method = "REML"
                    )
summary(mod_lnRR_mr)
# Multivariate Meta-Analysis Model (k = 318; method: REML)
# 
#    logLik   Deviance        AIC        BIC       AICc   
# -106.8735   213.7470   221.7470   236.7699   221.8756   
# 
# Variance Components:
# 
#             estim    sqrt  nlvls  fixed          factor 
# sigma^2.1  0.0614  0.2477     36     no        study_ID 
# sigma^2.2  0.0832  0.2885    318     no  effect_size_ID 
# 
# Test for Residual Heterogeneity:
# QE(df = 316) = 13892.9906, p-val < .0001
# 
# Test of Moderators (coefficient 2):
# F(df1 = 1, df2 = 316) = 5.7785, p-val = 0.0168
# 
# Model Results:
# 
#                           estimate      se     tval   df    pval    ci.lb    ci.ub      
# intrcpt                    -0.2663  0.0552  -4.8269  316  <.0001  -0.3749  -0.1578  *** 
# Org.fertilizer.typeplant   -0.1407  0.0585  -2.4038  316  0.0168  -0.2558  -0.0255    * 

r2_ml(mod_lnRR_mr)
# R2_marginal R2_conditional 
#   0.0323952      0.4429578 

The intercept represents the estimated mean lnRR for studies using animal-based fertilisers. The coefficient for plant-based fertilisers quantifies the difference in lnRR between plant- and animal-based fertilisers.

The intercept represents the estimated mean lnRR for studies using animal-based fertilisers, and the coefficient for plant-based fertilisers quantifies the difference in lnRR between plant- and animal-based fertilisers. The estimated coefficient for plant-based fertilisers was negative and statistically significant, indicating that yields under plant-based organic fertilisers were, on average, lower than those under animal-based organic fertilisers. The moderator explained approximately 3.2% of the total variance in effect sizes (\(R^2_{\text{marginal}}\) = 0.032) - it indicates that differences between fertiliser types account for a small portion of the overall heterogeneity in effect sizes.

\(I^2\) vs \(R^2\)

\(I^2\) and \(R^2\) are both measures used in meta-analysis, but they quantify different aspects of variability in effect sizes.

\(I^2\) measures the proportion of total variability in effect sizes that is due to heterogeneity among studies rather than sampling error. It provides insight into how much of the observed variation can be attributed to real differences between studies.

\(R^2\), on the other hand, quantifies the proportion of variance in effect sizes that is explained by the moderators included in a meta-regression model. It indicates how well the included predictors account for the variability in effect sizes. \(R^2_{\text{marginal}}\) measures the power of your chosen moderators to explain heterogeneity, while \(R^2_{\text{conditional}}\) measures the total variance accounted for by the model.

brms

As we mentioned in the main text, we can use the brms package to fit the same multilevel meta-analysis model. There are different ways to run these models, but in this example, we use the gr() function to set the sampling variance-covariance matrix. This approach is less common, but it can be much faster than usual (see more detailed information: link).

If you are not familiar with brms, we highly recommend checking out the instruction of brms first.

In the model below, sampling error is explicitly represented using a group-level term with a known variance–covariance structure. The term (1 | gr(effect_size_ID, cov = vcv)) specifies a random effect whose covariance matrix is fixed to the diagonal matrix of sampling variances for lnRR. This formulation is mathematically equivalent to supplying a sampling variance matrix in metafor::rma.mv(), but implemented here using brms.

Because sampling variances are assumed to be known in meta-analysis, the standard deviation of this group-level effect is fixed to 1. As a result, the magnitude of sampling error is fully determined by the supplied variance–covariance matrix and is not estimated from the data.

vcv <- diag(dat$var.lnRR)

rownames(vcv) <- colnames(vcv) <- dat$effect_size_ID

form_ma1 <- bf(lnRR ~ 1 +
                 (1|study_ID) + # this is u_l (the between-study effect)
                 (1|gr(effect_size_ID, cov = vcv)) # this is m (sampling error)
               )

# generate default priors
prior_ma1 <- default_prior(form_ma1, 
                           data=dat, 
                           data2=list(vcv=vcv),
                           family=gaussian())
prior_ma1$prior[3] = "constant(1)" # meta-analysis assumes sampling variance is known so fixing this to 1
prior_ma1

# fitting model
fit_ma1 <- brm(
  formula = form_ma1,
  data = dat,
  data2 = list(vcv=vcv),
  chains = 2,
  cores = 2,
  iter = 6000,
  warmup = 3000,
  prior = prior_ma1,
  control = list(adapt_delta=0.95, max_treedepth=15)
  )
summary(fit_ma1)
#  Family: gaussian 
#   Links: mu = identity 
# Formula: lnRR ~ 1 + (1 | study_ID) + (1 | gr(effect_size_ID, cov = vcv)) 
#    Data: dat (Number of observations: 318) 
#   Draws: 2 chains, each with iter = 6000; warmup = 3000; thin = 1;
#          total post-warmup draws = 6000
# 
# Multilevel Hyperparameters:
# ~effect_size_ID (Number of levels: 318) 
#               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
# sd(Intercept)     1.00      0.00     1.00     1.00   NA       NA       NA
# 
# ~study_ID (Number of levels: 36) 
#               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
# sd(Intercept)     0.28      0.04     0.21     0.38 1.00     1839     3325
# 
# Regression Coefficients:
#           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
# Intercept    -0.32      0.06    -0.43    -0.21 1.00     1842     2814
# 
# Further Distributional Parameters:
#       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
# sigma     0.29      0.01     0.26     0.32 1.00     7189     4944
# 
# Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
# and Tail_ESS are effective sample size measures, and Rhat is the potential
# scale reduction factor on split chains (at convergence, Rhat = 1).

The estimated overall mean and between-study standard deviation is comparable to the study-level variance component obtained from the metafor model. In contrast, the residual parameter sigma represents unexplained heterogeneity at the effect-size level beyond sampling error, corresponding to the effect-size-level variance component in rma.mv.

One important point concerns how random effects are reported. In brms, random effect estimates are presented as standard deviations (SD), whereas metafor reports variances . Therefore, to compare estimates between the two packages, you need to square the standard deviation from brms to obtain the variance.

Let’s now fit a meta-regression model including fertiliser type as a moderator - you can find the outputs are very similar to those from metafor.

vcv <- diag(dat$var.lnRR)
rownames(vcv) <- colnames(vcv) <- dat$effect_size_ID

form_mr1 <- bf(lnRR ~ 1 + Org.fertilizer.type +
                  (1|study_ID) + # this is u_l (the between-study effect)
                  (1|gr(effect_size_ID, cov = vcv)) # this is m (sampling error)
               )

# generate default priors
prior_mr1 <- default_prior(form_mr1,
                           data=dat, 
                           data2=list(vcv=vcv),
                           family=gaussian())
prior_mr1$prior[5] = "constant(1)" # meta-analysis assumes

# fitting model
fit_mr1 <- brm(
  formula = form_mr1,
  data = dat,
  data2 = list(vcv=vcv),
  chains = 2,
  cores = 2,
  iter = 6000,
  warmup = 3000,
  prior = prior_mr1,
  control = list(adapt_delta=0.95, max_treedepth=15)
  )
summary(fit_mr1)
#  Family: gaussian 
#   Links: mu = identity 
# Formula: lnRR ~ 1 + Org.fertilizer.type + (1 | study_ID) + (1 | gr(effect_size_ID, cov = vcv)) 
#    Data: dat (Number of observations: 318) 
#   Draws: 2 chains, each with iter = 6000; warmup = 3000; thin = 1;
#          total post-warmup draws = 6000
# 
# Multilevel Hyperparameters:
# ~effect_size_ID (Number of levels: 318) 
#               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
# sd(Intercept)     1.00      0.00     1.00     1.00   NA       NA       NA
# 
# ~study_ID (Number of levels: 36) 
#               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
# sd(Intercept)     0.26      0.04     0.19     0.35 1.00     2384     4173
# 
# Regression Coefficients:
#                          Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
# Intercept                   -0.27      0.06    -0.38    -0.15 1.00     2787     3503
# Org.fertilizer.typeplant    -0.14      0.06    -0.26    -0.02 1.00     4423     4398
# 
# Further Distributional Parameters:
#       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
# sigma     0.29      0.01     0.26     0.32 1.00     6530     4799
# 
# Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
# and Tail_ESS are effective sample size measures, and Rhat is the potential
# scale reduction factor on split chains (at convergence, Rhat = 1).

You may wonder why we specify ~1 + Org.fertilizer.type rather than simply ~ Org.fertilizer.type. In fact, there is no difference between these two specifications. In R model formulas, an intercept term is included by default, so both formulas fit a model with an intercept.

We write 1 explicitly here for clarity. Doing so makes it clear that the model includes an overall intercept, which represents the baseline mean effect size for the reference level of the moderator. This is particularly helpful in a meta-analytic context, where the intercept has a clear interpretation as the overall mean effect. For this reason, we explicitly include 1 in all model formulas throughout this tutorial. This convention makes the role of the intercept transparent and helps avoid confusion when comparing different model specifications.

glmmTMB

glmmTMB is the popular package for fitting generalised linear mixed models. Unlike dedicated meta-analytic packages such as metafor, it does not automatically account for known sampling variances of effect sizes - so when we implement a meta-analysis in glmmTMB, we need explicitly specify the sampling error structure of the effect sizes.

In practice, this require three steps:

1. Constructing a variance-covariance matrix (VCV) that contains the known sampling variances of each effect size

vcv <- diag(dat$var.lnRR)

row.names(vcv) <- colnames(vcv) <- dat$effect_size_ID

Here, we create a variance–covariance matrix whose diagonal elements are the known sampling variances of the effect sizes (var.lnRR). Off-diagonal elements are set to zero, reflecting the assumption that effect sizes are independent. This mirrors the standard assumption made in most meta-analyses.

2. Treating each effect size as a unique random-effect level

dat$effect_size_ID <- factor(dat$effect_size_ID)

In glmmTMB, grouping variables for random effects must be factors. By converting effect_size_ID to a factor, we ensure that each effect size is treated as a distinct level in the random-effects structure. This step is required for the model to correctly associate each observation with its corresponding sampling variance.

3. Constraining the corresponding random-effect variance structure to be equal to the known VCV

g <- rep(1, nrow(dat))

The equalto() term in glmmTMB requires a grouping factor, even when all observations belong to a single group. The dummy variable g serves this purpose. All observations are assigned to the same group, while the variance structure across effect sizes is fully defined by the supplied VCV matrix.

Conceptually, this step allows us to include sampling error as a random effect with a fixed, known variance structure.

You can also refer to this supplementary material for more details: https://coraliewilliams.github.io/equalto_sim_study/webpage.html.

glmm_lnRR_ma <- glmmTMB(lnRR ~ 1 +
                        (1 | study_ID) + 
                         equalto(0 + effect_size_ID|g, vcv), #  incorporates the known sampling variances of each effect size
                         data = dat,
                        REML = TRUE
                        )
summary(glmm_lnRR_ma)
#  Family: gaussian  ( identity )
# Formula:          lnRR ~ 1 + (1 | study_ID) + equalto(0 + effect_size_ID | g, vcv)
# Data: dat
# 
#       AIC       BIC    logLik -2*log(L)  df.resid
#     231.4     242.7    -112.7     225.4       315
# 
# Random effects:
# 
# Conditional model:
#  Groups   Name             Variance  Std.Dev. Corr
#  study_ID (Intercept)      0.0735694 0.27124
#  g        effect_size_ID33 0.0047213 0.06871
#           effect_size_ID34 0.0426425 0.20650  0.00
#           effect_size_ID35 0.0698313 0.26426  0.00 0.00
#           effect_size_ID36 0.0798567 0.28259  0.00 0.00 0.00
#           effect_size_ID37 0.0837435 0.28938  0.00 0.00 0.00 0.00
#           effect_size_ID39 0.0095870 0.09791  0.00 0.00 0.00 0.00 0.00
#           effect_size_ID40 0.0091520 0.09567  0.00 0.00 0.00 0.00 0.00 0.00
#           effect_size_ID42 0.0363323 0.19061  0.00 0.00 0.00 0.00 0.00 0.00 0.00
#           effect_size_ID44 0.0331452 0.18206  0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
#           effect_size_ID45 0.0012338 0.03513  0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
#           ...              ...       ...      ...  ...  ...  ...  ...  ...  ...  ...  ...  ...
#  Residual                  0.0830247 0.28814
# Number of obs: 318, groups:  study_ID, 36; g, 1
# 
# Dispersion estimate for gaussian family (sigma^2): 0.083
# 
# Conditional model:
#             Estimate Std. Error z value Pr(>|z|)
# (Intercept) -0.32039    0.05358   -5.98 2.23e-09 ***

In short, you should interpret the output as follows:

  • Overall effect estimate: (Intercept)

  • Between-study variance estimate: study_ID (Intercept) variance

  • Effect-size–level (sampling) variance: Residual / sigma^2

All estimate is identical to metafor output, and also close to brms one.

NOTE:
The grouping factor g is a technical dummy variable required by glmmTMB when using equalto() to impose a known variance-covariance structure. All observations belong to a single level of g, and the variance structure across effect sizes is entirely defined by the supplied VCV matrix.

As a result, the multiple variance components listed under g in the model output correspond to the known sampling variances of individual effect sizes and are not estimated parameters. These components should not be interpreted. Only the study-level variance and the residual variance are of substantive interest.

When you want to run the meta-regression, you will simply add the moderators in the model:

glmm_lnRR_mr <- glmmTMB(lnRR ~ 1 + Org.fertilizer.type + (1 | study_ID) + 
                      equalto(0 + effect_size_ID|g, vcv),
                  data = dat,
                  REML = TRUE
                  )
summary(glmm_lnRR_mr)
#  Family: gaussian  ( identity )
# Formula:          lnRR ~ 1 + Org.fertilizer.type + (1 | study_ID) + equalto(0 +  
#     effect_size_ID | g, vcv)
# Data: dat
# 
#       AIC       BIC    logLik -2*log(L)  df.resid 
#     231.9     246.9    -111.9     223.9       314 
# 
# Random effects:
# 
# Conditional model:
#  Groups   Name             Variance  Std.Dev. Corr                                            
#  study_ID (Intercept)      0.0613541 0.24770                                                  
#  g        effect_size_ID33 0.0047213 0.06871                                                  
#           effect_size_ID34 0.0426425 0.20650  0.00                                            
#           effect_size_ID35 0.0698313 0.26426  0.00 0.00                                       
#           effect_size_ID36 0.0798567 0.28259  0.00 0.00 0.00                                  
#           effect_size_ID37 0.0837435 0.28938  0.00 0.00 0.00 0.00                             
#           effect_size_ID39 0.0095870 0.09791  0.00 0.00 0.00 0.00 0.00                        
#           effect_size_ID40 0.0091520 0.09567  0.00 0.00 0.00 0.00 0.00 0.00                   
#           effect_size_ID42 0.0363323 0.19061  0.00 0.00 0.00 0.00 0.00 0.00 0.00              
#           effect_size_ID44 0.0331452 0.18206  0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00         
#           effect_size_ID45 0.0012338 0.03513  0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00    
#           ...              ...       ...      ...  ...  ...  ...  ...  ...  ...  ...  ...  ...
#  Residual                  0.0832439 0.28852                                                  
# Number of obs: 318, groups:  study_ID, 36; g, 1
# 
# Dispersion estimate for gaussian family (sigma^2): 0.0832 
# 
# Conditional model:
#                          Estimate Std. Error z value Pr(>|z|)    
# (Intercept)              -0.26633    0.05522  -4.823 1.41e-06 ***
# Org.fertilizer.typeplant -0.14069    0.06008  -2.342   0.0192 * 

Great - the estimates match those from metafor and are also very close to the results from brms :)

Model 2: Meta-analysis of variability

In this section we analyse variability effects using lnCVR and lnVR.

  • lnCVR quantifies the log ratio of coefficients of variation (CVs), that is, changes in relative variability after accounting for mean differences.

  • lnVR quantifies the log ratio of standard deviations (SDs), that is, changes in absolute variability between groups.

metafor

lnCVR

mod_lnCVR <- rma.mv(lnCVR,
                   var.lnCVR,
                   random = list(~ 1 | study_ID,
                                 ~ 1 | effect_size_ID),
                   data = dat,
                   test = "t",
                   method = "REML"
                   )
summary(mod_lnCVR)
# Multivariate Meta-Analysis Model (k = 318; method: REML)
# 
#    logLik   Deviance        AIC        BIC       AICc   
# -439.8288   879.6576   885.6576   896.9343   885.7342   
# 
# Variance Components:
# 
#             estim    sqrt  nlvls  fixed          factor 
# sigma^2.1  0.0602  0.2455     36     no        study_ID 
# sigma^2.2  0.2976  0.5455    318     no  effect_size_ID 
# 
# Test for Heterogeneity:
# Q(df = 317) = 803.2915, p-val < .0001
# 
# Model Results:
# 
# estimate      se    tval   df    pval   ci.lb   ci.ub      
#   0.3430  0.0709  4.8354  317  <.0001  0.2034  0.4825  *** 
# 
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

i2_ml(mod_lnCVR)
# I2_Total       I2_study_ID I2_effect_size_ID 
# 56.525070          9.517436         47.007634

The intercept from the multilevel meta-analysis represents the overall mean lnCVR across all comparisons. A positive estimate indicates that the treatment group has a larger CVthan the control group, meaning greater relative variability (lower stability) after accounting for potential mean differences.

Here, the estimated lnCVR is \(\beta_0\) = 0.34 (95% CI= [0.20, 0.48]), suggesting that organic farming is associated with approximately 41% higher relative variability in crop yields compared to conventional farming (= exp(0.34)−1).

\(I^2\) values indicate that a substantial portion of the total heterogeneity (\(I^2_{\text{Total}}\) = 56.53%) arises from differences between effect sizes within studies (\(I^2_{\text{within}}\) = 47.01%), while a smaller share is due to variability between studies (\(I^2_{\text{between}}\) = 9.52%).

lnVR

mod_lnVR <- rma.mv(lnVR,
                   var.lnVR,
                   random = list(~ 1 | study_ID,
                                 ~ 1 | effect_size_ID),
                   data = dat,
                   test = "t",
                   method = "REML"
                   )
summary(mod_lnVR)
# Multivariate Meta-Analysis Model (k = 318; method: REML)
# 
#    logLik   Deviance        AIC        BIC       AICc   
# -423.3326   846.6651   852.6651   863.9418   852.7418   
# 
# Variance Components:
# 
#             estim    sqrt  nlvls  fixed          factor 
# sigma^2.1  0.0640  0.2530     36     no        study_ID 
# sigma^2.2  0.1886  0.4342    318     no  effect_size_ID 
# 
# Test for Heterogeneity:
# Q(df = 317) = 658.8664, p-val < .0001
# 
# Model Results:
# 
# estimate      se    tval   df    pval    ci.lb   ci.ub    
#   0.0198  0.0676  0.2924  317  0.7702  -0.1133  0.1528    
# 
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

i2_ml(mod_lnVR)
# I2_Total       I2_study_ID I2_effect_size_ID 
# 50.62070          12.82949          37.79121 

The overall lnVR estimate is close to zero and the 95% CI overlaps zero (\(\beta_0\) = 0.02, 95% CI = [-0.11, 0.15]), indicating little evidence for a systematic difference in SD between the two groups. On the original scale, this corresponds to an SD ratio close to one, suggesting that the magnitude of within-group variation is broadly similar across conditions.

Despite the absence of a clear average lnVR effect, heterogeneity is substantial. The test for heterogeneity is highly significant, and approximately half of the observed variability in lnVR cannot be attributed to sampling error alone (\(I^2 \approx\) 50.6%). Only a smaller proportion of this heterogeneity is attributable to between-study differences (\(I^2 \approx\) 12.8%), whereas a larger share arises at the effect-size level within studies (\(I^2 \approx\) 37.8%). This pattern indicates considerable context dependence or within-study variation in absolute variability effects, even though the average effect across all comparisons is close to zero.

When considered alongside lnCVR, a clearer picture emerges. While lnVR shows little evidence for a consistent difference in absolute variability, lnCVR indicates a positive average effect, suggesting higher relative variability (CV) in the treatment group. This divergence arises because lnCVR explicitly standardises variability by the mean, whereas lnVR does not. When group means differ, absolute variability may remain similar (lnVR \(\approx\) 0) while relative variability increases (lnCVR > 0).

In practical terms, this pattern implies that outcomes may scale with the mean. Although the absolute spread of values is comparable across groups, the spread is larger relative to the expected outcome under the treatment condition. Consequently, differences in variability are more apparent on the mean-standardised scale than on the absolute scale.

lnCVR

mod_lnCVR_mr <- rma.mv(lnCVR,
                    var.lnCVR,
                    mods = ~ Org.fertilizer.type,
                    random = list(~ 1 | study_ID,
                                  ~ 1 | effect_size_ID),
                    data = dat,
                    test = "t",
                    method = "REML"
                    )
summary(mod_lnCVR_mr)
# Multivariate Meta-Analysis Model (k = 318; method: REML)
# 
#    logLik   Deviance        AIC        BIC       AICc   
# -438.7285   877.4570   885.4570   900.4800   885.5856   
# 
# Variance Components:
# 
#             estim    sqrt  nlvls  fixed          factor 
# sigma^2.1  0.0602  0.2454     36     no        study_ID 
# sigma^2.2  0.3001  0.5478    318     no  effect_size_ID 
# 
# Test for Residual Heterogeneity:
# QE(df = 316) = 803.2270, p-val < .0001
# 
# Test of Moderators (coefficient 2):
# F(df1 = 1, df2 = 316) = 0.3026, p-val = 0.5826
# 
# Model Results:
# 
#                           estimate      se    tval   df    pval    ci.lb   ci.ub      
# intrcpt                     0.3150  0.0874  3.6020  316  0.0004   0.1429  0.4870  *** 
# Org.fertilizer.typeplant    0.0644  0.1171  0.5501  316  0.5826  -0.1660  0.2949      

r2_ml(mod_lnCVR_mr)
   # R2_marginal R2_conditional 
   #  0.00281066     0.16943708 

lnVR

mod_lnVR_mr <- rma.mv(lnVR,
                      var.lnVR,
                      mods = ~ Org.fertilizer.type,
                      random = list(~ 1 | study_ID,
                                    ~ 1 | effect_size_ID),
                      data = dat,
                      test = "t",
                      method = "REML"
                      )
summary(mod_lnVR_mr)
# Multivariate Meta-Analysis Model (k = 318; method: REML)
# 
#    logLik   Deviance        AIC        BIC       AICc   
# -421.3055   842.6109   850.6109   865.6339   850.7396   
# 
# Variance Components:
# 
#             estim    sqrt  nlvls  fixed          factor 
# sigma^2.1  0.0560  0.2367     36     no        study_ID 
# sigma^2.2  0.1896  0.4355    318     no  effect_size_ID 
# 
# Test for Residual Heterogeneity:
# QE(df = 316) = 651.1955, p-val < .0001
# 
# Test of Moderators (coefficient 2):
# F(df1 = 1, df2 = 316) = 2.4353, p-val = 0.1196
# 
# Model Results:
# 
#                           estimate      se     tval   df    pval    ci.lb   ci.ub    
# intrcpt                     0.0876  0.0791   1.1075  316  0.2689  -0.0680  0.2433    
# Org.fertilizer.typeplant   -0.1644  0.1054  -1.5605  316  0.1196  -0.3717  0.0429    

r2_ml(mod_lnVR_mr)
# R2_marginal R2_conditional 
#  0.02621268     0.24827290 

Including fertiliser type as a moderator provided little evidence that it systematically explains variability effects on either scale. Here, the moderator coefficient represents the difference between plant-based and animal-based fertilisers (plant − animal).

For lnCVR, the estimated moderator effect was small and highly uncertain (\(\beta_{\text{plant-animal}}\) = 0.064, 95% CI = [-0.17, 0.29]), indicating that relative variability does not differ consistently between fertiliser categories.

Similarly, for lnVR, there was no clear support for fertiliser-type differences in absolute variability (\(\beta_{\text{plant-animal}}\) = −0.164, 95% CI = [-0.37, 0.04]). Although the point estimate suggested a tendency toward lower lnVR for plant-based fertilisers relative to animal-based fertilisers, this effect was highly uncertain and not statistically supported.

In both models, substantial residual heterogeneity remained after accounting for fertiliser type, with large effect-size–level variance components (lnCVR: \(\sigma^2\) = 0.30; lnVR: \(\sigma^2\) = 0.19). This indicates that neither relative nor absolute variability effects are well captured by this single moderator.

brms

When fitted with equivalent model structures, the metafor and brms implementations yield nearly identical estimates for both the overall effects and variance components.

lnCVR

vcv <- diag(dat$var.lnCVR)

rownames(vcv) <- colnames(vcv) <- dat$effect_size_ID

form_ma2 <- bf(lnCVR ~ 1 +
                 (1|study_ID) + # this is u_l (the between-study effect)
                 (1|gr(effect_size_ID, cov = vcv)) # this is m (sampling error)
               )

# generate default priors
prior_ma2 <- default_prior(form_ma2, 
                           data=dat, 
                           data2=list(vcv=vcv),
                           family=gaussian())
prior_ma2$prior[3] = "constant(1)" # meta-analysis assumes sampling variance is known so fixing this to 1

# fitting model
fit_ma2 <- brm(
  formula = form_ma2,
  data = dat,
  data2 = list(vcv=vcv),
  chains = 2,
  cores = 2,
  iter = 6000,
  warmup = 3000,
  prior = prior_ma2,
  control = list(adapt_delta=0.95, max_treedepth=15)
  )

summary(fit_ma2)
# Family: gaussian 
#   Links: mu = identity 
# Formula: lnCVR ~ 1 + (1 | study_ID) + (1 | gr(effect_size_ID, cov = vcv)) 
#    Data: dat (Number of observations: 318) 
#   Draws: 2 chains, each with iter = 6000; warmup = 3000; thin = 1;
#          total post-warmup draws = 6000
# 
# Multilevel Hyperparameters:
# ~effect_size_ID (Number of levels: 318) 
#               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
# sd(Intercept)     1.00      0.00     1.00     1.00   NA       NA       NA
# 
# ~study_ID (Number of levels: 36) 
#               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
# sd(Intercept)     0.27      0.08     0.12     0.45 1.00     1404     2356
# 
# Regression Coefficients:
#           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
# Intercept     0.34      0.07     0.20     0.49 1.00     3888     4078
# 
# Further Distributional Parameters:
#       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
# sigma     0.55      0.05     0.45     0.66 1.00      896     1820

lnVR

vcv <- diag(dat$var.lnVR)

rownames(vcv) <- colnames(vcv) <- dat$effect_size_ID

form_ma3 <- bf(lnVR ~ 1 +
                 (1|study_ID) + # this is u_l (the between-study effect)
                 (1|gr(effect_size_ID, cov = vcv)) # this is m (sampling error)
               )

# generate default priors
prior_ma3 <- default_prior(form_ma3, 
                           data=dat, 
                           data2=list(vcv=vcv),
                           family=gaussian())
prior_ma3$prior[3] = "constant(1)" # meta-analysis assumes sampling variance is known so fixing this to 1

# fitting model
fit_ma3 <- brm(
  formula = form_ma3,
  data = dat,
  data2 = list(vcv=vcv),
  chains = 2,
  cores = 2,
  iter = 6000,
  warmup = 3000,
  prior = prior_ma3,
  control = list(adapt_delta=0.95, max_treedepth=15)
  )
summary(fit_ma3)
#  Family: gaussian 
#   Links: mu = identity 
# Formula: lnVR ~ 1 + (1 | study_ID) + (1 | gr(effect_size_ID, cov = vcv)) 
#    Data: dat (Number of observations: 318) 
#   Draws: 2 chains, each with iter = 6000; warmup = 3000; thin = 1;
#          total post-warmup draws = 6000
# 
# Multilevel Hyperparameters:
# ~effect_size_ID (Number of levels: 318) 
#               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
# sd(Intercept)     1.00      0.00     1.00     1.00   NA       NA       NA
# 
# ~study_ID (Number of levels: 36) 
#               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
# sd(Intercept)     0.27      0.07     0.15     0.42 1.00     1458     2196
# 
# Regression Coefficients:
#           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
# Intercept     0.02      0.07    -0.12     0.17 1.00     2816     3765
# 
# Further Distributional Parameters:
#       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
# sigma     0.44      0.06     0.32     0.56 1.01      469      856

When specified with equivalent fixed and random effects and the same sampling variance structure, the brms and metafor meta-regression models yield nearly identical results. Differences between outputs reflect estimation framework (Bayesian vs frequentist) rather than substantive differences in model inference.

lnCVR

vcv <- diag(dat$var.lnCVR)

rownames(vcv) <- colnames(vcv) <- dat$effect_size_ID

form_mr2 <- bf(lnCVR ~ 1 + Org.fertilizer.type +
                 (1|study_ID) + # this is u_l (the between-study effect)
                 (1|gr(effect_size_ID, cov = vcv)) # this is m (sampling error)
               )

# generate default priors
prior_mr2 <- default_prior(form_mr2,
                           data=dat, 
                           data2=list(vcv=vcv),
                           family=gaussian())
prior_mr2$prior[5] = "constant(1)" # meta-analysis assumes

# fitting model
fit_mr2 <- brm(
  formula = form_mr2,
  data = dat,
  data2 = list(vcv=vcv),
  chains = 2,
  cores = 2,
  iter = 6000,
  warmup = 3000,
  prior = prior_mr2,
  control = list(adapt_delta=0.95, max_treedepth=15)
  )
summary(fit_mr2)
#  Family: gaussian 
#   Links: mu = identity 
# Formula: lnCVR ~ 1 + Org.fertilizer.type + (1 | study_ID) + (1 | gr(effect_size_ID, cov = vcv)) 
#    Data: dat (Number of observations: 318) 
#   Draws: 2 chains, each with iter = 6000; warmup = 3000; thin = 1;
#          total post-warmup draws = 6000
# 
# Multilevel Hyperparameters:
# ~effect_size_ID (Number of levels: 318) 
#               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
# sd(Intercept)     1.00      0.00     1.00     1.00   NA       NA       NA
# 
# ~study_ID (Number of levels: 36) 
#               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
# sd(Intercept)     0.27      0.08     0.13     0.43 1.00     1515     2605
# 
# Regression Coefficients:
#                          Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
# Intercept                    0.32      0.09     0.14     0.50 1.00     4506     4337
# Org.fertilizer.typeplant     0.06      0.12    -0.18     0.30 1.00     3709     4471
# 
# Further Distributional Parameters:
#       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
# sigma     0.55      0.05     0.45     0.66 1.00      987     1853

lnVR

vcv <- diag(dat$var.lnVR)

rownames(vcv) <- colnames(vcv) <- dat$effect_size_ID

form_mr3 <- bf(lnVR ~ 1 + Org.fertilizer.type +
                 (1|study_ID) + # this is u_l (the between-study effect)
                 (1|gr(effect_size_ID, cov = vcv)) # this is m (sampling error)
               )

# generate default priors
prior_mr3 <- default_prior(form_mr3,
                           data=dat, 
                           data2=list(vcv=vcv),
                           family=gaussian())
prior_mr3$prior[5] = "constant(1)" # meta-analysis assumes

# fitting model
fit_mr3 <- brm(
  formula = form_mr3,
  data = dat,
  data2 = list(vcv=vcv),
  chains = 2,
  cores = 2,
  iter = 6000,
  warmup = 3000,
  prior = prior_mr3,
  control = list(adapt_delta=0.95, max_treedepth=15)
  )
summary(fit_mr3)
# Family: gaussian 
#   Links: mu = identity 
# Formula: lnVR ~ 1 + Org.fertilizer.type + (1 | study_ID) + (1 | gr(effect_size_ID, cov = vcv)) 
#    Data: dat (Number of observations: 318) 
#   Draws: 2 chains, each with iter = 6000; warmup = 3000; thin = 1;
#          total post-warmup draws = 6000
# 
# Multilevel Hyperparameters:
# ~effect_size_ID (Number of levels: 318) 
#               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
# sd(Intercept)     1.00      0.00     1.00     1.00   NA       NA       NA
# 
# ~study_ID (Number of levels: 36) 
#               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
# sd(Intercept)     0.24      0.07     0.11     0.39 1.00     1038     1472
# 
# Regression Coefficients:
#                          Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
# Intercept                    0.09      0.08    -0.06     0.25 1.00     3110     4251
# Org.fertilizer.typeplant    -0.17      0.11    -0.38     0.04 1.00     2781     3712
# 
# Further Distributional Parameters:
#       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
# sigma     0.44      0.06     0.32     0.56 1.00      443      659

glmmTMB

We can also fit the same meta-analytic and meta-regression models using glmmTMB. You can find that the outputs are very similar to those from metafor and brms.

First, we need to prepare the variance-covariance matrix for the sampling errors as shown below…

# lnCVR
vcv <- diag(dat$var.lnCVR)
row.names(vcv) <- colnames(vcv) <- dat$effect_size_ID
dat$effect_size_ID <- factor(dat$effect_size_ID)
g <- rep(1, nrow(dat))

# lnVR
vcv <- diag(dat$var.lnVR)
row.names(vcv) <- colnames(vcv) <- dat$effect_size_ID
dat$effect_size_ID <- factor(dat$effect_size_ID)
g <- rep(1, nrow(dat))
glmm_lnCVR <- glmmTMB(lnCVR ~ 1 + (1 | study_ID) + 
                        equalto(0 + effect_size_ID|g, vcv),
                      data = dat,
                      REML = TRUE
                      )
summary(glmm_lnCVR)
#  Family: gaussian  ( identity )
# Formula:          lnCVR ~ 1 + (1 | study_ID) + equalto(0 + effect_size_ID | g,      vcv)
# Data: dat
# 
#       AIC       BIC    logLik -2*log(L)  df.resid 
#     892.6     903.8    -443.3     886.6       315 
# 
# Random effects:
# 
# Conditional model:
#  Groups   Name             Variance Std.Dev. Corr                                            
#  study_ID (Intercept)      0.06173  0.2485                                                   
#  g        effect_size_ID33 0.11111  0.3333                                                   
#           effect_size_ID34 0.25000  0.5000   0.00                                            
#           effect_size_ID35 0.25000  0.5000   0.00 0.00                                       
#           effect_size_ID36 0.25000  0.5000   0.00 0.00 0.00                                  
#           effect_size_ID37 0.25000  0.5000   0.00 0.00 0.00 0.00                             
#           effect_size_ID39 0.11111  0.3333   0.00 0.00 0.00 0.00 0.00                        
#           effect_size_ID40 0.11111  0.3333   0.00 0.00 0.00 0.00 0.00 0.00                   
#           effect_size_ID42 0.25000  0.5000   0.00 0.00 0.00 0.00 0.00 0.00 0.00              
#           effect_size_ID44 0.25000  0.5000   0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00         
#           effect_size_ID45 0.25000  0.5000   0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00    
#           ...              ...      ...      ...  ...  ...  ...  ...  ...  ...  ...  ...  ...
#  Residual                  0.31923  0.5650                                                   
# Number of obs: 318, groups:  study_ID, 36; g, 1
# 
# Dispersion estimate for gaussian family (sigma^2): 0.319 
# 
# Conditional model:
#             Estimate Std. Error z value Pr(>|z|)    
# (Intercept)  0.34418    0.07145   4.817 1.46e-06 ***
glmm_lnVR <- glmmTMB(lnVR ~ 1 + (1 | study_ID) + 
                        equalto(0 + effect_size_ID|g, vcv),
                      data = dat,
                      REML = TRUE
                     )
summary(glmm_lnVR)
#  Family: gaussian  ( identity )
# Formula:          lnVR ~ 1 + (1 | study_ID) + equalto(0 + effect_size_ID | g, vcv)
# Data: dat
# 
#       AIC       BIC    logLik -2*log(L)  df.resid 
#     858.4     869.7    -426.2     852.4       315 
# 
# Random effects:
# 
# Conditional model:
#  Groups   Name             Variance Std.Dev. Corr                                            
#  study_ID (Intercept)      0.06401  0.2530                                                   
#  g        effect_size_ID33 0.11111  0.3333                                                   
#           effect_size_ID34 0.25000  0.5000   0.00                                            
#           effect_size_ID35 0.25000  0.5000   0.00 0.00                                       
#           effect_size_ID36 0.25000  0.5000   0.00 0.00 0.00                                  
#           effect_size_ID37 0.25000  0.5000   0.00 0.00 0.00 0.00                             
#           effect_size_ID39 0.11111  0.3333   0.00 0.00 0.00 0.00 0.00                        
#           effect_size_ID40 0.11111  0.3333   0.00 0.00 0.00 0.00 0.00 0.00                   
#           effect_size_ID42 0.25000  0.5000   0.00 0.00 0.00 0.00 0.00 0.00 0.00              
#           effect_size_ID44 0.25000  0.5000   0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00         
#           effect_size_ID45 0.25000  0.5000   0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00    
#           ...              ...      ...      ...  ...  ...  ...  ...  ...  ...  ...  ...  ...
#  Residual                  0.18856  0.4342                                                   
# Number of obs: 318, groups:  study_ID, 36; g, 1
# 
# Dispersion estimate for gaussian family (sigma^2): 0.189 
# 
# Conditional model:
#             Estimate Std. Error z value Pr(>|z|)
# (Intercept)  0.01977    0.06787   0.291    0.771
glmm_lnCVR_mr <- glmmTMB(lnCVR ~ 1
                        + Org.fertilizer.type
                        + (1 | study_ID) + 
                        equalto(0 + effect_size_ID|g, vcv),
                      data = dat,
                      REML = TRUE
                      )
summary(glmm_lnCVR_mr)
#  Family: gaussian  ( identity )
# Formula:          lnCVR ~ 1 + Org.fertilizer.type + (1 | study_ID) + equalto(0 +      effect_size_ID | g, vcv)
# Data: dat
# 
#       AIC       BIC    logLik -2*log(L)  df.resid 
#     896.6     911.7    -444.3     888.6       314 
# 
# Random effects:
# 
# Conditional model:
#  Groups   Name             Variance Std.Dev. Corr                                            
#  study_ID (Intercept)      0.06125  0.2475                                                   
#  g        effect_size_ID33 0.11111  0.3333                                                   
#           effect_size_ID34 0.25000  0.5000   0.00                                            
#           effect_size_ID35 0.25000  0.5000   0.00 0.00                                       
#           effect_size_ID36 0.25000  0.5000   0.00 0.00 0.00                                  
#           effect_size_ID37 0.25000  0.5000   0.00 0.00 0.00 0.00                             
#           effect_size_ID39 0.11111  0.3333   0.00 0.00 0.00 0.00 0.00                        
#           effect_size_ID40 0.11111  0.3333   0.00 0.00 0.00 0.00 0.00 0.00                   
#           effect_size_ID42 0.25000  0.5000   0.00 0.00 0.00 0.00 0.00 0.00 0.00              
#           effect_size_ID44 0.25000  0.5000   0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00         
#           effect_size_ID45 0.25000  0.5000   0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00    
#           ...              ...      ...      ...  ...  ...  ...  ...  ...  ...  ...  ...  ...
#  Residual                  0.32171  0.5672                                                   
# Number of obs: 318, groups:  study_ID, 36; g, 1
# 
# Dispersion estimate for gaussian family (sigma^2): 0.322 
# 
# Conditional model:
#                          Estimate Std. Error z value Pr(>|z|)    
# (Intercept)               0.31394    0.08819   3.560 0.000371 ***
# Org.fertilizer.typeplant  0.07016    0.11845   0.592 0.553664    
glmm_lnVR_mr <- glmmTMB(lnVR ~ 1
                        + Org.fertilizer.type
                        + (1 | study_ID) +
                        equalto(0 + effect_size_ID|g, vcv),
                      data = dat,
                      REML = TRUE
                      )
summary(glmm_lnVR_mr)
#  Family: gaussian  ( identity )
# Formula:          lnVR ~ 1 + Org.fertilizer.type + (1 | study_ID) + equalto(0 +      effect_size_ID | g, vcv)
# Data: dat
# 
#       AIC       BIC    logLik -2*log(L)  df.resid 
#     860.7     875.8    -426.4     852.7       314 
# 
# Random effects:
# 
# Conditional model:
#  Groups   Name             Variance Std.Dev. Corr                                            
#  study_ID (Intercept)      0.05602  0.2367                                                   
#  g        effect_size_ID33 0.11111  0.3333                                                   
#           effect_size_ID34 0.25000  0.5000   0.00                                            
#           effect_size_ID35 0.25000  0.5000   0.00 0.00                                       
#           effect_size_ID36 0.25000  0.5000   0.00 0.00 0.00                                  
#           effect_size_ID37 0.25000  0.5000   0.00 0.00 0.00 0.00                             
#           effect_size_ID39 0.11111  0.3333   0.00 0.00 0.00 0.00 0.00                        
#           effect_size_ID40 0.11111  0.3333   0.00 0.00 0.00 0.00 0.00 0.00                   
#           effect_size_ID42 0.25000  0.5000   0.00 0.00 0.00 0.00 0.00 0.00 0.00              
#           effect_size_ID44 0.25000  0.5000   0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00         
#           effect_size_ID45 0.25000  0.5000   0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00    
#           ...              ...      ...      ...  ...  ...  ...  ...  ...  ...  ...  ...  ...
#  Residual                  0.18965  0.4355                                                   
# Number of obs: 318, groups:  study_ID, 36; g, 1
# 
# Dispersion estimate for gaussian family (sigma^2): 0.19 
# 
# Conditional model:
#                          Estimate Std. Error z value Pr(>|z|)
# (Intercept)               0.08762    0.07921   1.106    0.269
# Org.fertilizer.typeplant -0.16443    0.10620  -1.548    0.122

Models 3 and 4: Location-scale meta-regression models

We now extend the above meta-analyses to location–scale meta-regression models, which allow us to examine whether both the mean effect size and the variability of effect sizes depend on a moderator variable. Here, we use Org.fertilizer.type (plant-based vs animal-based organic fertiliser) as a moderator.

  • The location component answers the question:
    Does the mean effect size (i.e. lnRR, lnCVR, lnVR) differ between plant-based and animal-based organic fertilisers?

  • The scale component answers the question:
    Does the residual heterogeneity (i.e. variability of effect sizes after accounting for random effects) differ between plant-based and animal-based organic fertilisers?

metafor

At present, metafor does not support fully multilevel location–scale meta-analytic models. However, it is possible to fit a simpler location–scale model using the scale argument in rma(). This approach allows the residual heterogeneity to vary as a function of moderators, but does not accommodate multilevel random-effects structures.

Below, we demonstrate two types of models:

  1. Simple location–scale meta-regression

  2. Multilevel heteroscedastic meta-regression

Note

It is best to treat the metafor models presented above as approximations to full multilevel location-scale meta-analytic models. Although these models are extremely useful and sometimes sufficient in practice, ecological datasets often show more complex hierarchical structures, such as multiple effect sizes per study, evolutionary independence, or other sources of non-independence. To explicitly consider such complexity, more flexible modelling frameworks are required. In the following two sections, we introduce location-scale meta-analytic models fitted using brms and glmmTMB. These approach arrow the specification of rich multilevel structures, and in particular, brms enables the fitting of fully Bayesian location-scale models with complex random effects and covariance structures (though if we want to run this type of model, we need a lot of data…)

The formulation using rma() provides a straightforward way to test whether variability in effect sizes depends on the moderator, but it treats all effect sizes as independent and does not include a multilevel structure. scale argument allows us to specify a model for the residual variance, enabling us to test whether the variability of effect sizes differs across levels of the moderator.

lnRR

mod_lnRR_ls <- rma(yi = lnRR,
                   vi = var.lnRR,
                   mods = ~ Org.fertilizer.type,
                   scale = ~ Org.fertilizer.type,
                   data = dat,
                   test = "t",
                   method = "REML"
                   )
summary(mod_lnRR_ls)
# Location-Scale Model (k = 318; tau^2 estimator: REML)
# 
#    logLik   deviance        AIC        BIC       AICc   
# -161.2053   322.4106   330.4106   345.4336   330.5393   
# 
# Test for Residual Heterogeneity:
# QE(df = 316) = 13892.9906, p-val < .0001
# 
# Test of Location Coefficients (coefficient 2):
# F(df1 = 1, df2 = 316) = 67.0786, p-val < .0001
# 
# Model Results (Location):
# 
#                           estimate      se     tval   df    pval    ci.lb    ci.ub      
# intrcpt                    -0.1796  0.0290  -6.1935  316  <.0001  -0.2366  -0.1225  *** 
# Org.fertilizer.typeplant   -0.3666  0.0448  -8.1902  316  <.0001  -0.4547  -0.2785  *** 
# 
# Test of Scale Coefficients (coefficient 2):
# F(df1 = 1, df2 = 316) = 19.9223, p-val < .0001
# 
# Model Results (Scale):
# 
#                           estimate      se      tval   df    pval    ci.lb    ci.ub      
# intrcpt                    -2.4057  0.1423  -16.9060  316  <.0001  -2.6856  -2.1257  *** 
# Org.fertilizer.typeplant    0.8071  0.1808    4.4634  316  <.0001   0.4513   1.1629  *** 

lnCVR

mod_lnCVR_ls <- rma(yi = lnCVR,
                   vi = var.lnCVR,
                   mods = ~ Org.fertilizer.type,
                   scale = ~ Org.fertilizer.type,
                   data = dat,
                   test = "t",
                   method = "REML"
                   )
summary(mod_lnCVR_ls)
# Location-Scale Model (k = 318; tau^2 estimator: REML)
# 
#    logLik   deviance        AIC        BIC       AICc   
# -439.5156   879.0311   887.0311   902.0541   887.1597   
# 
# Test for Residual Heterogeneity:
# QE(df = 316) = 803.2270, p-val < .0001
# 
# Test of Location Coefficients (coefficient 2):
# F(df1 = 1, df2 = 316) = 1.0750, p-val = 0.3006
# 
# Model Results (Location):
# 
#                           estimate      se    tval   df    pval    ci.lb   ci.ub      
# intrcpt                     0.2401  0.0620  3.8719  316  0.0001   0.1181  0.3621  *** 
# Org.fertilizer.typeplant    0.1026  0.0990  1.0368  316  0.3006  -0.0921  0.2974      
# 
# Test of Scale Coefficients (coefficient 2):
# F(df1 = 1, df2 = 316) = 10.1534, p-val = 0.0016
# 
# Model Results (Scale):
# 
#                           estimate      se     tval   df    pval    ci.lb    ci.ub      
# intrcpt                    -1.4908  0.2473  -6.0286  316  <.0001  -1.9773  -1.0043  *** 
# Org.fertilizer.typeplant    1.0399  0.3263   3.1864  316  0.0016   0.3978   1.6820   ** 

lnVR

mod_lnVR_ls <- rma(yi = lnVR,
                   vi = var.lnVR,
                   mods = ~ Org.fertilizer.type,
                   scale = ~ Org.fertilizer.type,
                   data = dat,
                   test = "t",
                   method = "REML"
                   )
summary(mod_lnVR_ls)
# Location-Scale Model (k = 318; tau^2 estimator: REML)
# 
#    logLik   deviance        AIC        BIC       AICc   
# -415.2505   830.5010   838.5010   853.5240   838.6296   
# 
# Test for Residual Heterogeneity:
# QE(df = 316) = 651.1955, p-val < .0001
# 
# Test of Location Coefficients (coefficient 2):
# F(df1 = 1, df2 = 316) = 4.7895, p-val = 0.0294
# 
# Model Results (Location):
# 
#                           estimate      se     tval   df    pval    ci.lb    ci.ub    
# intrcpt                     0.0417  0.0457   0.9106  316  0.3632  -0.0483   0.1317    
# Org.fertilizer.typeplant   -0.1906  0.0871  -2.1885  316  0.0294  -0.3620  -0.0193  * 
# 
# Test of Scale Coefficients (coefficient 2):
# F(df1 = 1, df2 = 316) = 11.9043, p-val = 0.0006
# 
# Model Results (Scale):
# 
#                           estimate      se     tval   df    pval    ci.lb    ci.ub      
# intrcpt                    -2.8652  0.6270  -4.5699  316  <.0001  -4.0987  -1.6316  *** 
# Org.fertilizer.typeplant    2.3044  0.6679   3.4503  316  0.0006   0.9903   3.6185  *** 

The models suggest that plant-based fertilisers are associated with both lower mean effects and substantially greater variability across studies. In particular, the scale component indicates that effect sizes are markedly more dispersed for plant-based fertilisers across all three variability metrics. However, these models treat all effect sizes as independent and do not account for the multilevel structure of the data.

Heteroscedastic meta-regression models allow residual variance to differ across levels of a moderator while incorporating a multilevel random-effects structure. By relaxing the assumption of constant heterogeneity, these models permit different levels of a moderator to be associated with different amounts of residual variability. Although this approach does not constitute a full multilevel location–scale model, it provides a useful intermediate step toward more flexible variance modelling. In particular, it enables the examination of moderator effects on both the mean effect size and the residual heterogeneity while accounting for the multilevel structure of the data.

In the present heteroscedastic meta-regression, the random-effects structure ~ Org.fertilizer.type | effect_size_ID allows residual heterogeneity at the effect-size level to differ across fertiliser types. By specifying struct = "DIAG", the model estimates fertiliser-type-specific variances without assuming covariance between categories, thereby relaxing the assumption of homogeneous residual variability while retaining a parsimonious specification.

lnRR

mod_lnRR_mr2 <- rma.mv(lnRR,
                       var.lnRR,
                       mods = ~ Org.fertilizer.type,
                       random = list(~ 1 | study_ID,
                                     ~ Org.fertilizer.type | effect_size_ID),
                       data = dat,
                       test = "t",
                       method = "REML",
                       struct =  'DIAG'
                       )
summary(mod_lnRR_mr2)
# Multivariate Meta-Analysis Model (k = 318; method: REML)
# 
#    logLik   Deviance        AIC        BIC       AICc   
# -105.4125   210.8250   220.8250   239.6037   221.0185   
# 
# Variance Components:
# 
#             estim    sqrt  nlvls  fixed    factor 
# sigma^2    0.0608  0.2466     36     no  study_ID 
# 
# outer factor: effect_size_ID      (nlvls = 318)
# inner factor: Org.fertilizer.type (nlvls = 2)
# 
#             estim    sqrt  k.lvl  fixed   level 
# tau^2.1    0.0666  0.2581    134     no  animal 
# tau^2.2    0.0937  0.3061    184     no   plant 
# 
# Test for Residual Heterogeneity:
# QE(df = 316) = 13892.9906, p-val < .0001
# 
# Test of Moderators (coefficient 2):
# F(df1 = 1, df2 = 316) = 5.9453, p-val = 0.0153
# 
# Model Results:
# 
#                           estimate      se     tval   df    pval    ci.lb    ci.ub      
# intrcpt                    -0.2679  0.0537  -4.9852  316  <.0001  -0.3736  -0.1622  *** 
# Org.fertilizer.typeplant   -0.1399  0.0574  -2.4383  316  0.0153  -0.2528  -0.0270    * 

lnCVR

mod_lnCVR_mr2 <- rma.mv(lnCVR,
                       var.lnCVR,
                       mods = ~ Org.fertilizer.type,
                       random = list(~ 1 | study_ID,
                                     ~ Org.fertilizer.type | effect_size_ID),
                       data = dat,
                       test = "t",
                       method = "REML",
                       struct =  'DIAG'
                       )
summary(mod_lnCVR_mr2)
# Multivariate Meta-Analysis Model (k = 318; method: REML)
# 
#    logLik   Deviance        AIC        BIC       AICc   
# -434.6673   869.3347   879.3347   898.1134   879.5282   
# 
# Variance Components:
# 
#             estim    sqrt  nlvls  fixed    factor 
# sigma^2    0.0578  0.2405     36     no  study_ID 
# 
# outer factor: effect_size_ID      (nlvls = 318)
# inner factor: Org.fertilizer.type (nlvls = 2)
# 
#             estim    sqrt  k.lvl  fixed   level 
# tau^2.1    0.1915  0.4376    134     no  animal 
# tau^2.2    0.5419  0.7361    184     no   plant 
# 
# Test for Residual Heterogeneity:
# QE(df = 316) = 803.2270, p-val < .0001
# 
# Test of Moderators (coefficient 2):
# F(df1 = 1, df2 = 316) = 0.4957, p-val = 0.4819
# 
# Model Results:
# 
#                           estimate      se    tval   df    pval    ci.lb   ci.ub      
# intrcpt                     0.3021  0.0821  3.6816  316  0.0003   0.1406  0.4635  *** 
# Org.fertilizer.typeplant    0.0845  0.1200  0.7040  316  0.4819  -0.1517  0.3207      

lnVR

mod_lnVR_mr2 <- rma.mv(lnVR,
                       var.lnVR,
                       mods = ~ Org.fertilizer.type,
                       random = list(~ 1 | study_ID,
                                     ~ Org.fertilizer.type | effect_size_ID),
                       data = dat,
                       test = "t",
                       method = "REML",
                       struct =  'DIAG'
                       )
summary(mod_lnVR_mr2)
# Multivariate Meta-Analysis Model (k = 318; method: REML)
# 
#    logLik   Deviance        AIC        BIC       AICc   
# -409.8610   819.7221   829.7221   848.5008   829.9156   
# 
# Variance Components:
# 
#             estim    sqrt  nlvls  fixed    factor 
# sigma^2    0.0588  0.2425     36     no  study_ID 
# 
# outer factor: effect_size_ID      (nlvls = 318)
# inner factor: Org.fertilizer.type (nlvls = 2)
# 
#             estim    sqrt  k.lvl  fixed   level 
# tau^2.1    0.0291  0.1707    134     no  animal 
# tau^2.2    0.5092  0.7136    184     no   plant 
# 
# Test for Residual Heterogeneity:
# QE(df = 316) = 651.1955, p-val < .0001
# 
# Test of Moderators (coefficient 2):
# F(df1 = 1, df2 = 316) = 1.9172, p-val = 0.1671
# 
# Model Results:
# 
#                           estimate      se     tval   df    pval    ci.lb   ci.ub    
# intrcpt                     0.0666  0.0713   0.9329  316  0.3516  -0.0738  0.2069    
# Org.fertilizer.typeplant   -0.1526  0.1102  -1.3846  316  0.1671  -0.3695  0.0643    

In the heteroscedastic meta-regression for lnRR, residual heterogeneity at the effect-size level differed between fertiliser types, with larger heterogeneity for plant-based fertilisers (\(\tau^2\) = 0.094) than for animal-based fertilisers (\(\tau^2\) = 0.067). The location component indicated a significantly lower mean lnRR for plant-based fertilisers relative to animal-based fertilisers (estimate = −0.14, 95% CI = [-0.25, -0.03]). Thus, even after accounting for multilevel structure, plant-based fertilisers were associated with both lower mean effects and moderately greater residual heterogeneity.

For lnCVR, allowing residual heterogeneity to vary across fertiliser types revealed substantial differences in heterogeneity magnitude. Effect-size–level heterogeneity was markedly larger for plant-based fertilisers (\(\tau^2\) = 0.54) than for animal-based fertilisers (\(\tau^2\) = 0.19). However, the location effect of fertiliser type on mean lnCVR remained small and statistically unsupported (estimate = 0.08, 95% CI = [-0.15, 0.32]), indicating no consistent difference in average relative variability between fertiliser categories.

A similar pattern was observed for lnVR. Residual heterogeneity was considerably larger for plant-based fertilisers (\(\tau^2\) = 0.51) than for animal-based fertilisers (\(\tau^2\) = 0.03), suggesting pronounced differences in the dispersion of absolute variability across fertiliser types. In contrast, the location effect on mean lnVR was not statistically supported (estimate = −0.15, 95% CI = [-0.37, 0.06]), indicating weak evidence for systematic differences in average absolute variability.

brms

brms provides a flexible framework for fitting location–scale (LS) models, in which the mean structure of the response variable (the location component) and the residual variability around that mean (the scale component) are modelled simultaneously. This allows us to examine not only how predictors affect average effect sizes, but also how they influence heterogeneity.

In this tutorial, we focus primarily on a location–scale meta-regression model with fixed effects in the scale component. This model represents the most stable and interpretable formulation and serves as the main working example throughout the brms section.

For completeness, we also briefly introduce two more complex extensions:

  1. Location–scale meta-regression with fixed effects in the scale component This is the primary model used in this tutorial and the only one for which we provide full interpretation of results.

  2. Location–scale meta-regression with random effects in the scale component This model allows heterogeneity to vary across higher-level grouping factors (for example, studies), but may be difficult to fit reliably in practice.

  3. Location–scale meta-regression with correlated random effects in the location and scale components This formulation allows the mean and heterogeneity to be evolutionarily or structurally linked, but often suffers from convergence issues in applied meta-analytic datasets.

Models 2 and 3 are included to illustrate the range of model structures that are theoretically possible in brms. We provide one concrete example of each, but we do not interpret their results in detail, as these models frequently exhibit convergence problems and require a lot of data. If you are interested in more detail, please check Nakagawa et al. (2025a, b).

It is worth noting that, in the models presented here, we specify the scale component using the sigma ~ syntax in brms. This choice reflects the use of a Gaussian response distribution, for which sigma represents the residual standard deviation. For other response distributions, the scale component may be parametrised differently in brms (for example, using shape or phi - see here), and the corresponding model formula should be adapted accordingly.

Practical guidance

Models with random effects in the scale component, and especially models with correlated random effects across location and scale, require:

  • Sufficient numbers of studies

  • Adequate replication within studies

  • Careful prior specification

  • Thorough convergence diagnostics

For these reasons, we treat Models 2 and 3 as examples. We recommend that readers begin with a simpler location–scale model that uses fixed effects in the scale part, and only add complexity if the data clearly support it, and also, you have enough data.

In this model, we include fixed effects in both location and scale components, but no random effects in the scale component (of course, we can include any random effects on location component).

Between-study heterogeneity in mean effect sizes is captured through a random intercept for study_ID in the location component, and residual heterogeneity is allowed to vary systematically with predictors via fixed effects in the scale component. No random variation in residual variance across studies is assumed.

lnRR

form_ls_ma1 <- bf(lnRR ~ 1 + Org.fertilizer.type +
                    (1|study_ID) + # this is u_l (the between-study effect)
                    (1|gr(effect_size_ID, cov = vcv)), # this is m (sampling error)
                  sigma ~ 1 +  Org.fertilizer.type # fixed effects only in scale component
                  )

# generate default priors
prior_ls_ma1 <- default_prior(form_ls_ma1, 
                              data=dat, 
                              data2=list(vcv=vcv),
                              family=gaussian())

prior_ls_ma1$prior[5] = "constant(1)" # meta-analysis assumes

# fitting model
fit_ls_ma1 <- brm(
  formula = form_ls_ma1,
  data = dat,
  data2 = list(vcv=vcv),
  chains = 2,
  cores = 2,
  iter = 6000,
  warmup = 3000,
  prior = prior_ls_ma1,
  control = list(adapt_delta=0.95, max_treedepth=15)
  )
summary(fit_ls_ma1)
# Family: gaussian 
#   Links: mu = identity; sigma = log 
# Formula: lnRR ~ 1 + Org.fertilizer.type + (1 | study_ID) + (1 | gr(effect_size_ID, cov = vcv)) 
#          sigma ~ 1 + Org.fertilizer.type
#    Data: dat (Number of observations: 318) 
#   Draws: 2 chains, each with iter = 6000; warmup = 3000; thin = 1;
#          total post-warmup draws = 6000
# 
# Multilevel Hyperparameters:
# ~effect_size_ID (Number of levels: 318) 
#               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
# sd(Intercept)     0.45      0.27     0.02     1.02 1.00      865     2021
# 
# ~study_ID (Number of levels: 36) 
#               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
# sd(Intercept)     0.26      0.04     0.19     0.35 1.00     2342     3505
# 
# Regression Coefficients:
#                                Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
# Intercept                         -0.26      0.06    -0.38    -0.16 1.00     2437     3607
# sigma_Intercept                   -1.30      0.07    -1.45    -1.15 1.00     3773     3686
# Org.fertilizer.typeplant          -0.15      0.06    -0.27    -0.03 1.00     5617     4995
# sigma_Org.fertilizer.typeplant     0.17      0.09    -0.01     0.35 1.00     9838     4318

For lnRR, the location component indicated a negative average effect under animal-based fertilisers (\(\beta^{(l)}_{{\text{animal}}\) = -0.26, 95% credible interval CrI = [−0.38, −0.16]), implying that organic yields are, on average, lower than conventional yields under animal-based fertilisation. The contrast coefficient for fertiliser type was also negative (\(\beta^{(l)}_{\text{animal-plant}}\) = -0.15, 95% CrI [−0.27, −0.03]), indicating an additional yield penalty under plant-based fertilisers relative to animal-based fertilisers. On the response-ratio scale, these estimates correspond approximately to organic yields of (-0.26) \(\approx\) 0.77 (about 23% lower than conventional) under animal-based fertilisers, versus (-0.26 - 0.15) \(\approx\) 0.66 (about 34% lower) under plant-based fertilisers.

The scale component, modelled on the log-\(\sigma\) scale, suggested a tendency towards greater residual heterogeneity under plant-based fertilisers (\(\beta^{(s)}_{\text{animal-plant}}\) = 0.17, 95% CrI = [−0.01, 0.35]). This corresponds to an approximate 19% increase in residual heterogeneity for plant-based systems ((0.17)-1 \(\approx\) 0.19), although the credible interval still overlapped zero, indicating moderate uncertainty in the magnitude of this effect.

lnCVR

form_ls_ma2 <- bf(lnCVR ~ 1 + Org.fertilizer.type +
                    (1|study_ID) + # this is u_l (the between-study effect)
                    (1|gr(effect_size_ID, cov = vcv)), # this is m (sampling error)
                  sigma ~ 1 +  Org.fertilizer.type
                  )

# generate default priors
prior_ls_ma2 <- default_prior(form_ls_ma2, 
                              data=dat, 
                              data2=list(vcv=vcv),
                              family=gaussian())

prior_ls_ma2$prior[5] = "constant(1)" # meta-analysis assumes

# fitting model
fit_ls_ma2 <- brm(
  formula = form_ls_ma2,
  data = dat,
  data2 = list(vcv=vcv),
  chains = 2,
  cores = 2,
  iter = 6000,
  warmup = 3000,
  prior = prior_ls_ma2,
  control = list(adapt_delta=0.95, max_treedepth=15)
  )
summary(fit_ls_ma2)
# Family: gaussian 
#   Links: mu = identity; sigma = log 
# Formula: lnCVR ~ 1 + Org.fertilizer.type + (1 | study_ID) + (1 | gr(effect_size_ID, cov = vcv)) 
#          sigma ~ 1 + Org.fertilizer.type
#    Data: dat (Number of observations: 318) 
#   Draws: 2 chains, each with iter = 6000; warmup = 3000; thin = 1;
#          total post-warmup draws = 6000
# 
# Multilevel Hyperparameters:
# ~effect_size_ID (Number of levels: 318) 
#               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
# sd(Intercept)     1.00      0.00     1.00     1.00   NA       NA       NA
# 
# ~study_ID (Number of levels: 36) 
#               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
# sd(Intercept)     0.26      0.08     0.12     0.43 1.00     1516     2572
# 
# Regression Coefficients:
#                                Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
# Intercept                          0.30      0.09     0.13     0.48 1.00     4592     4409
# sigma_Intercept                   -0.83      0.13    -1.09    -0.58 1.00     1518     2937
# Org.fertilizer.typeplant           0.08      0.13    -0.17     0.32 1.00     4734     4478
# sigma_Org.fertilizer.typeplant     0.51      0.18     0.15     0.87 1.00     1459     2679
# 
# Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
# and Tail_ESS are effective sample size measures, and Rhat is the potential
# scale reduction factor on split chains (at convergence, Rhat = 1).

For lnCVR, the location component indicated that organic systems tend to exhibit higher mean-corrected variability under animal-based fertilisers (\(\beta^{(l)}_{\text{animal}}\) = 0.30, 95% CrI = [0.13, 0.48]), corresponding to an average CV ratio of (0.30) \(\approx\) 1.35 - roughly 35% higher variability in organic relative to conventional systems.

The contrast coefficient for fertiliser type in the location component was small and highly uncertain (\(\beta^{(l)}_{\text{animal-plant}}\) = 0.08, 95% CrI [−0.17, 0.32]), suggesting that fertiliser category does not systematically shift the average stability difference between organic and conventional systems.

In contrast, the scale component revealed clear evidence that fertiliser type strongly governs the predictability of lnCVR effects. Residual heterogeneity was substantially larger under plant-based fertilisers \((\beta^{(s)}_{\text{animal-plant}}\) = 0.51, 95% CrI [0.15, 0.87]), corresponding to a multiplicative increase in residual variability of approximately 68% ((0.51) \(\approx\) 1.68). Under this model, the implied residual standard deviation for lnCVR under plant-based fertilisers is \(\sigma \approx \exp(-0.83 + 0.51) \approx 0.73\).

Thus, while the expected magnitude of the stability gap (lnCVR) is broadly similar across fertiliser categories, the evidence is far less transferable under plant-based fertilisers. Some studies show near parity in stability between organic and conventional systems, whereas others show much larger differences, indicating strong context dependence.

The models presented here extend the basic location–scale meta-regression by allowing study-level heterogeneity to operate not only on the mean effect size, but also on the residual variance.

These models are included for completeness and to illustrate the flexibility of brms. In practice, they are often difficult to fit reliably and should be treated as advanced extensions rather than default choices. Accordingly, we provide example code but do not interpret the results in detail. If you want to see the examples, please refer to this online tutorial.

2. with random effects in scale component

In the model with the random effects in the scale component, we can include random effects in the scale component, allowing residual heterogeneity to vary across studies. If we want to specify such a model, we can do so as follows:

sigma ~ 1 + Org.fertilizer.type + (1 | study_ID)

This specification allows the residual variance to vary across studies, in addition to allowing study-level variation in the mean effect size via the location component. Conceptually, the location component models how the mean effect size varies with the moderator and across studies. The scale component models how the amount of unexplained variability differs across studies and fertiliser types.

Random effects in the location and scale components are assumed to be independent in this formulation.

3. Correlated random effects in location and scale components

The third model further extends the previous formulation by allowing correlated random effects between the location and scale components. For example, if we want to consider potential correlations between study-level deviations in mean effect sizes and deviations in residual variance, we can specify the model as follows:

(1 | p | study_ID)

This specification allows the model to estimate the correlation between study-level deviations in the mean effect size and deviations in residual variance.

For example, a negative correlation would indicate that studies with larger mean effects tend to show lower residual variability, whereas a positive correlation would suggest the opposite. Such relationships are conceptually appealing and biologically interpretable in some contexts.

However, this model is computationally demanding and highly sensitive to prior specification and data structure. In applied meta-analyses, it frequently exhibits convergence problems or weakly identified correlation parameters.

glmmTMB

In contrast to metafor, glmmTMB allows us to include random effects in both location and scale components - we demonstrate how to fit a location-scale meta-regression model using glmmTMB.

But it should be mentioned that glmmTMB does not support correlated random effects between location and scale components, so the third model presented above cannot be implemented using glmmTMB.

Again, the same as brms, we focus on a location–scale meta-regression model with fixed effects in the scale component.

1. Only fixed effects in scale component

In glmmTMB, the location–scale structure is implemented by separating the model into a conditional component and a dispersion component. The conditional model specifies the mean structure of the response (location), while the dispersion model specifies how the residual variance changes as a function of predictors (scale). Note that the dispersion model in glmmTMB does not describe an additional random effect, but rather models heteroscedasticity in the residual variance on the log scale.

As shown below, the location–scale models with glmmTMB produce results that match those from brms completely.

lnRR

glmm_lnRR1 <- glmmTMB(lnRR ~ 1 + Org.fertilizer.type + (1 | study_ID) + 
                      equalto(0 + effect_size_ID|g, vcv),
                  dispformula = ~ 1 + Org.fertilizer.type, # scale component
                  data = dat,
                  REML = TRUE)
summary(glmm_lnRR1)
#  Family: gaussian  ( identity )
# Formula:          lnRR ~ 1 + Org.fertilizer.type + (1 | study_ID) + equalto(0 +      effect_size_ID | g, vcv)
# Dispersion:            ~1 + Org.fertilizer.type
# Data: dat
# 
#       AIC       BIC    logLik -2*log(L)  df.resid 
#     230.9     249.7    -110.5     220.9       313 
# 
# Random effects:
# 
# Conditional model:
#  Groups   Name             Variance  Std.Dev. Corr                                            
#  study_ID (Intercept)      0.0608129 0.24660                                                  
#  g        effect_size_ID33 0.0047213 0.06871                                                  
#           effect_size_ID34 0.0426425 0.20650  0.00                                            
#           effect_size_ID35 0.0698313 0.26426  0.00 0.00                                       
#           effect_size_ID36 0.0798567 0.28259  0.00 0.00 0.00                                  
#           effect_size_ID37 0.0837435 0.28938  0.00 0.00 0.00 0.00                             
#           effect_size_ID39 0.0095870 0.09791  0.00 0.00 0.00 0.00 0.00                        
#           effect_size_ID40 0.0091520 0.09567  0.00 0.00 0.00 0.00 0.00 0.00                   
#           effect_size_ID42 0.0363323 0.19061  0.00 0.00 0.00 0.00 0.00 0.00 0.00              
#           effect_size_ID44 0.0331452 0.18206  0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00         
#           effect_size_ID45 0.0012338 0.03513  0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00    
#           ...              ...       ...      ...  ...  ...  ...  ...  ...  ...  ...  ...  ...
#  Residual                         NA      NA                                                  
# Number of obs: 318, groups:  study_ID, 36; g, 1
# 
# Conditional model:
#                          Estimate Std. Error z value Pr(>|z|)    
# (Intercept)              -0.26791    0.05379  -4.981 6.33e-07 ***
# Org.fertilizer.typeplant -0.13991    0.05897  -2.373   0.0177 *  
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 
# Dispersion model:
#                          Estimate Std. Error z value Pr(>|z|)    
# (Intercept)              -1.35445    0.07922 -17.097   <2e-16 ***
# Org.fertilizer.typeplant  0.17055    0.09928   1.718   0.0858 . 

lnCVR

glmm_lnCVR1 <- glmmTMB(lnCVR ~ 1 + Org.fertilizer.type + (1 | study_ID) + 
                        equalto(0 + effect_size_ID|g, vcv),
                      dispformula = ~ 1 + Org.fertilizer.type,
                      data = dat,
                      REML = TRUE)
summary(glmm_lnCVR1)
#  Family: gaussian  ( identity )
# Formula:          lnCVR ~ 1 + Org.fertilizer.type + (1 | study_ID) + equalto(0 +      effect_size_ID | g, vcv)
# Dispersion:             ~1 + Org.fertilizer.type
# Data: dat
# 
#       AIC       BIC    logLik -2*log(L)  df.resid 
#     889.4     908.3    -439.7     879.4       313 
# 
# Random effects:
# 
# Conditional model:
#  Groups   Name             Variance Std.Dev. Corr                                            
#  study_ID (Intercept)      0.05782  0.2405                                                   
#  g        effect_size_ID33 0.11583  0.3403                                                   
#           effect_size_ID34 0.29264  0.5410   0.00                                            
#           effect_size_ID35 0.31983  0.5655   0.00 0.00                                       
#           effect_size_ID36 0.32986  0.5743   0.00 0.00 0.00                                  
#           effect_size_ID37 0.33374  0.5777   0.00 0.00 0.00 0.00                             
#           effect_size_ID39 0.12070  0.3474   0.00 0.00 0.00 0.00 0.00                        
#           effect_size_ID40 0.12026  0.3468   0.00 0.00 0.00 0.00 0.00 0.00                   
#           effect_size_ID42 0.28633  0.5351   0.00 0.00 0.00 0.00 0.00 0.00 0.00              
#           effect_size_ID44 0.28315  0.5321   0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00         
#           effect_size_ID45 0.25123  0.5012   0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00    
#           ...              ...      ...      ...  ...  ...  ...  ...  ...  ...  ...  ...  ...
#  Residual                       NA      NA                                                   
# Number of obs: 318, groups:  study_ID, 36; g, 1
# 
# Conditional model:
#                          Estimate Std. Error z value Pr(>|z|)    
# (Intercept)               0.30208    0.08297   3.641 0.000272 ***
# Org.fertilizer.typeplant  0.08451    0.12103   0.698 0.484986    
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 
# Dispersion model:
#                          Estimate Std. Error z value Pr(>|z|)    
# (Intercept)               -0.8265     0.1300  -6.356 2.07e-10 ***
# Org.fertilizer.typeplant   0.5201     0.1767   2.943  0.00325 ** 

2. Random effects in location and scale components

If you want to include random effects in the scale component, you just need to insert them to the dispformula part of the model.

glmm_lnRR2 <- glmmTMB(lnRR ~ 1 + Org.fertilizer.type + (1 | study_ID) + 
                        equalto(0 + effect_size_ID|g, vcv),
                      dispformula = ~ 1 + Org.fertilizer.type + (1 | study_ID),
                      data = dat,
                      REML = TRUE)

Gallary

Sometimes, visualising the outputs of meta-analysis models can be challenging, especially for brms models. Here, we show the scripts of visualisations for the meta-analysis and meta-regression models fitted using brms.

get_variables_dynamic <- function(model, patterns) {
  variables <- get_variables(model)
  idx <- unique(unlist(lapply(patterns, function(p) grep(p, variables))))
  if (length(idx) == 0) return(character(0))
  variables[idx]
}

# preparation
rename_vars <- function(variable) {
  # fixed effects
  variable <- gsub("b_Intercept", "b_l_intercept", variable)
  variable <- gsub("b_sigma_Intercept", "b_s_intercept", variable)  
  variable <- gsub("b_Org.fertilizer.typeplant", "b_l_contrast", variable)
  variable <- gsub("b_sigma_Org.fertilizer.typeplant", "b_s_contrast", variable)            
  # random effects
  variable <- gsub("sd_study_ID__Intercept", "sd_study_ID", variable)  
  variable <- gsub("sigma", "sd_effect_ID", variable)  
  return(variable)
}

# Function to visualise fixed effects
visualize_fixed_effects <- function(model) {
  fixed_effect_vars <- get_variables_dynamic(model, "^b_")
  if (length(fixed_effect_vars) == 0) {
    message("No fixed effects found")
    return(NULL)
  }
  
  tryCatch({
    fixed_effects_samples <- model %>%
      spread_draws(!!!syms(fixed_effect_vars)) %>%
      pivot_longer(cols = all_of(fixed_effect_vars), names_to = ".variable", values_to = ".value") %>%
      mutate(.variable = rename_vars(.variable))
    
    ggplot(fixed_effects_samples, aes(x = .value, y = .variable)) +
      stat_halfeye(
        normalize = "xy", 
        point_interval = "mean_qi", 
        fill = "lightcyan3", 
        color = "lightcyan4"
      ) +
      geom_vline(xintercept = 0, linetype = "dashed", color = "#005") +
      labs(y = "Fixed effects", x = "Posterior values") +
      theme_classic()
  }, error = function(e) {
    message("Error in visualize_fixed_effects: ", e$message)
    return(NULL)
  })
}

# Function to visualize random effects
visualize_random_effects <- function(model) {
  random_effect_vars <- get_variables_dynamic(model, c("^sd_study", "^sigma"))
  #random_effect_vars <- random_effect_vars[random_effect_vars != "sd_es_ID__Intercept"]
  if (length(random_effect_vars) == 0) {
    message("No random effects found")
    return(NULL)
  }
  
  tryCatch({
    random_effects_samples <- model %>%
      spread_draws(!!!syms(random_effect_vars)) %>%
      pivot_longer(cols = all_of(random_effect_vars), names_to = ".variable", values_to = ".value") %>%
      mutate(.variable = rename_vars(.variable)) #%>%
    #mutate(.value = .value)  # leave SD as it is
    
    ggplot(random_effects_samples, aes(x = .value, y = .variable)) +
      stat_halfeye(
        normalize = "xy", 
        point_interval = "mean_qi", 
        fill = "olivedrab3", 
        color = "olivedrab4"
      ) +
      geom_vline(xintercept = 0, linetype = "dashed", color = "#005") +
      labs(y = "Random effects (SD)", x = "Posterior values") +
      theme_classic()
  }, error = function(e) {
    message("Error in visualize_random_effects: ", e$message)
    return(NULL)
  })
}

# Function to visualize random effects
visualize_random_effects2 <- function(model) {
  random_effect_vars <- get_variables_dynamic(model, c("^sd_study"))
  #random_effect_vars <- random_effect_vars[random_effect_vars != "sd_es_ID__Intercept"]
  if (length(random_effect_vars) == 0) {
    message("No random effects found")
    return(NULL)
  }
  
  tryCatch({
    random_effects_samples <- model %>%
      spread_draws(!!!syms(random_effect_vars)) %>%
      pivot_longer(cols = all_of(random_effect_vars), names_to = ".variable", values_to = ".value") %>%
      mutate(.variable = rename_vars(.variable)) #%>%
    #mutate(.value = .value)  # leave SD as it is
    
    ggplot(random_effects_samples, aes(x = .value, y = .variable)) +
      stat_halfeye(
        normalize = "xy", 
        point_interval = "mean_qi", 
        fill = "olivedrab3", 
        color = "olivedrab4"
      ) +
      geom_vline(xintercept = 0, linetype = "dashed", color = "#005") +
      labs(y = "Random effects (SD)", x = "Posterior values") +
      theme_classic()
  }, error = function(e) {
    message("Error in visualize_random_effects: ", e$message)
    return(NULL)
  })
}

Software and package versions

R version 4.5.2 (2025-10-31)
Platform: aarch64-apple-darwin20
Running under: macOS Sequoia 15.7.3

Matrix products: default
BLAS:   /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib 
LAPACK: /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.1

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

time zone: America/Edmonton
tzcode source: internal

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] patchwork_1.3.2     here_1.0.2          orchaRd_2.1.3      
 [4] tidybayes_3.0.7     lubridate_1.9.4     forcats_1.0.1      
 [7] stringr_1.6.0       dplyr_1.1.4         purrr_1.2.0        
[10] readr_2.1.6         tidyr_1.3.2         tibble_3.3.0       
[13] ggplot2_4.0.1       tidyverse_2.0.0     metafor_4.8-0      
[16] numDeriv_2016.8-1.1 metadat_1.4-0       Matrix_1.7-4       
[19] brms_2.23.0         Rcpp_1.1.0         

loaded via a namespace (and not attached):
 [1] tidyselect_1.2.1      svUnit_1.0.8          vipor_0.4.7          
 [4] farver_2.1.2          loo_2.9.0             S7_0.2.1             
 [7] fastmap_1.2.0         latex2exp_0.9.6       tensorA_0.36.2.1     
[10] pacman_0.5.1          mathjaxr_2.0-0        digest_0.6.39        
[13] timechange_0.3.0      estimability_1.5.1    lifecycle_1.0.4      
[16] StanHeaders_2.32.10   magrittr_2.0.4        posterior_1.6.1      
[19] compiler_4.5.2        rlang_1.1.6           tools_4.5.2          
[22] igraph_2.2.1          yaml_2.3.12           phangorn_2.12.1      
[25] knitr_1.51            labeling_0.4.3        bridgesampling_1.2-1 
[28] htmlwidgets_1.6.4     pkgbuild_1.4.8        RColorBrewer_1.1-3   
[31] abind_1.4-8           withr_3.0.2           grid_4.5.2           
[34] stats4_4.5.2          colorspace_2.1-2      xtable_1.8-4         
[37] inline_0.3.21         emmeans_2.0.1         scales_1.4.0         
[40] cli_3.6.5             mvtnorm_1.3-3         rmarkdown_2.30       
[43] generics_0.1.4        otel_0.2.0            RcppParallel_5.1.11-1
[46] rstudioapi_0.17.1     tzdb_0.5.0            ape_5.8-1            
[49] ggbeeswarm_0.7.3      rstan_2.32.7          bayesplot_1.15.0     
[52] parallel_4.5.2        matrixStats_1.5.0     vctrs_0.6.5          
[55] jsonlite_2.0.0        hms_1.1.4             arrayhelpers_1.1-0   
[58] beeswarm_0.4.0        ggdist_3.3.3          glue_1.8.0           
[61] codetools_0.2-20      distributional_0.5.0  stringi_1.8.7        
[64] gtable_0.3.6          QuickJSR_1.8.1        quadprog_1.5-8       
[67] pillar_1.11.1         htmltools_0.5.9       Brobdingnag_1.2-9    
[70] R6_2.6.1              rprojroot_2.1.1       evaluate_1.0.5       
[73] lattice_0.22-7        backports_1.5.0       rstantools_2.5.0     
[76] fastmatch_1.1-6       gridExtra_2.3         coda_0.19-4.1        
[79] nlme_3.1-168          checkmate_2.3.3       xfun_0.55            
[82] pkgconfig_2.0.3      

References

  • Lauren C. Ponisio, Leithen K. M’Gonigle, Kevi C. Mace, Jenny Palomino, Perry de Valpine, Claire Kremen; Diversification practices reduce organic to conventional yield gap. Proc Biol Sci. 2015; 282 (1799): 20141396. https://doi.org/10.1098/rspb.2014.1396